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Measuring the Doppler broadening of the positron annihilation radiation or the angular correlation 
between the two annihilation gamma quanta reflects the momentum distribution of electrons seen 
by positrons in the material. Vacancy-type defects in solids localize positrons and the measured 
spectra are sensitive to the detailed chemical and geometric environments of the defects. However, 
the measured information is indirect and when using it in defect identification comparisons with 
theoretically predicted spectra is indispensable. In this article we present a computational scheme 
for calculating momentum distributions of electron-positron pairs annihilating in solids. Valence 
electron states and their interaction with ion cores are described using the all-electron projector 
augmented-wave method, and atomic orbitals are used to describe the core states. We apply our 
numerical scheme to selected systems and compare three different enhancement (electron-positron 
correlation) schemes previously used in the calculation of momentum distributions of annihilating 
electron-positron pairs within the density-functional theory. We show that the use of a state- 
dependent enhancement scheme leads to better results than a position-dependent enhancement 
factor in the case of ratios of Doppler spectra between different systems. Further, we demonstrate 
the applicability of our scheme for studying vacancy-type defects in metals and semiconductors. 
Especially we study the effect of forces due to a positron localized at a vacancy-type defect on the 
ionic relaxations. 



I. INTRODUCTION 

Positron annihilation spectroscopjii is an experimental 
method for studying electronic structures of materials. 
In comparison with other techniques to measure elec- 
tron momentum densities such as (e, 2e) spectroscopjiS 
or Compton scattering^ (to which we have already ap- 
plied the same all-electron method as used in this work^) 
positron annihilation is characterized with the strong sen- 
sitivity to the vacancy-type defects, which makes it a 
method widely suitable in materials science and materi- 
als technology studies. In the crystal lattice positrons 
get trapped at possibly existing vacancy-type defects. 
By measuring positron lifetimes and momentum distri- 
butions of annihilating electron-positron pairs (angular 
correlation or Doppler broadening measurements of an- 
nihilation radiation) one obtains information about the 
open volumes and the chemical environments of the de- 
fects. 

A successful use of positron annihilation measurements 
(in defect identification) calls for accompanying theoret- 
ical and computational work resulting in simulated an- 
nihilation characteristics to be compared with the mea- 
sured ones (for a review see Ref. i). In this paper 
we present a computational scheme based on the zero- 
positron-density {n^ — > 0) limit of the two-component 
density- functional theory^ (TCDFT). We describe the 
valence electron states in materials using the projector 
augmented-wave (PAW) methodj^ which we have already 
used to simulate the electron momentum distributions 
measured by Compton scattering.* In the PAW method 
the core states are treated in the frozen-core approxima- 
tion and described using atomic wave functions^ In this 



work, the positron state is solved in the real space us- 
ing a Rayleigh quotient multigrid (RQMG) solver.^ Our 
scheme gives good results when compared to experiments 
for delocalized positron states, for which the n+ — > 
limit of the TCDFT is exact, as well as for positrons 
localized at vacancy-type defects. 

The methods previously used in self-consistent cal- 
culations of electronic structures and wave functions 
for determination of momentum distributions of anni- 
hilating valence-electron-positron pairs include, for ex- 
ample, the pseudopotential methodfi2*ii*iSiiii4 the full- 
potential linearized-augmented-plane-wave (FLAPW) 
methodfi^iiSiii the hnear muffin-tin orbital (LMTO) 
methodic and localized basis set schemesJ^ Each of 
these methods has its own advantages and flaws. The 
FLAPW method, although being accurate, is computa- 
tionally heavy and has a basis set that makes the momen- 
tum density calculations technically complicated. The 
pseudopotential method is efficient and simple to use in 
the momentum distribution calculations. The drawback 
is that one completely loses the information on the high- 
momentum Fourier components of the valence wave func- 
tions because soft pseudo wave functions are used in the 
calculation. The LMTO method has presentation prob- 
lems in the interstitial regions, which renders it difficult 
to describe open structures like vacancy-type defects and 
systems with low symmetry with it. The PAW method, 
in contrast, is as efficient as the ultrasoft pseudopotential 
method. It can flexibly be used for the study of defects 
in solids including structural relaxation. The plane-wave 
representation of the pseudo wave functions makes things 
simple because plane- waves are eigenfunctions of the mo- 
mentum. Moreover, the calculation of the PAW momen- 
tum density^ is straightforward because the plane-waves 
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extend also to the augmentation regions as opposed to 
the (L)APW method. Ishibashi has aheady apphed the 
PAW method to the calculation of coincidence Doppler 
spectra for bulk materials, and Uenodo et al. have used 
the same code to study vacancy-type defects in SiGcj^ 
Previously, we have used the PAW scheme (without tak- 
ing into account the effect of the positron-induced forces) 
to study vacancy-dopant complexes in highly Sb doped 
Si,^^ monovacancy in Al,^^ and to show that one can ex- 
perimentally distinguish Ga vacancies (Vca) from Vca" 
On pairs in GaNpSS, 

Many works, in which the interpretation of the results 
is based on the comparison between measured and sim- 
ulated annihilation characteristics, have been published 
during the recent years. Therefore, a systematic study of 
the performance of different schemes and approximation 
is of utmost importance. Using the PAW method to de- 
scribe the valence electrons and atomic orbitals for core 
electrons we test three different schemes and approxima- 
tions to calculate momentum distributions of annihilat- 
ing electron-positron pairs. 

For the description of the many-body effects in the 
calculation of momentum distributions of annihilat- 
ing electron-positron pairs we choose finally the so- 
called state-dependent scheme^ and for annihilation rates 
within the state-dependent scheme the local-density ap- 
proximation (LDA) enhancement factor parametrized by 
Borohski and Nieminen>2i We show that use of the com- 
monly used position-dependent enhancement factor leads 
to unphysical oscillations at high momenta when one con- 
siders ratios of Doppler spectra between two different 
materials. In the same way we compare the Boronski- 
Nieminen LDA (BN-LDA) to the generalized-gradient 
approximation (GGA) by Barbiellini et aZ.^*^-^^ Our re- 
sults show that the BN-LDA describes the ratios more 
accurately compared with the experiment than the GGA. 
We also show that the ratios can be compared with 
the experiment reliably although the LDA enhancement 
overestimates the high-momentum region of the Doppler 
spectra arising from the annihilation with core electrons. 
There are no test systems solved theoretically {e.g. by 
the Quantum Monte Carlo method) exactly enough to 
compare with. 

When studying annihilation of positrons trapped at 
vacancies and comparing the results with experiments it 
is important to consider the effects of forces due to the 
localized positron on the ionic relaxation of the vacancy. 
The effects on calculated positron lifetimes and Doppler 
spectra are non-negligible. We study selected monova- 
cancies in metals and semiconductors by including also 
the effects of the forces due to the localized positron. 

The structure of the paper is as follows. In section UTI 
we describe the computational method used. Section UTTI 
presents results for bulk systems and section Hvl selected 
results for vacancies in metals and semiconductors. Fi- 
nally, we summarize our results and present conclusions 
in section Ivl 



II. THEORY AND COMPUTATIONAL 
METHODS 

A. Calculation of the positron states 

In our scheme wc first calculate the self-consistent elec- 
tronic structure of the system without the influence of 
the positron. Then we solve the positron state in the 
potential 

K^(r)=0(r)+Ko„(n_(r)), (1) 

where 0(r) is the Coulomb potential due to electrons and 
nuclei, ri_(r) the electron density, and T4orr("--(r)) is the 
—>■ limit of the electron-positron correlation poten- 
tial. Above, the self-interaction correction is explicitly 
made, i.e. since we are interested in the case of only one 
positron in the system, its self-direct Coulomb energy 
should exactly cancel its exchange-correlation energy. 
The calculation of the positron state is non-selfconsistent 
because the effective potential for the positron [Eq. ^] 
does not depend on the positron density. 

The scheme described above is for a delocalized 
positron the exact n+ ^ limit of the TCDFT but it 
has been proven appropriate in practice also in the case 
of localized positron states. In this connection it is of- 
ten referred to as the "conventional scheme"! (CONV). 
The approximation can be justified by considering the 
positron and its screening cloud as a neutral quasipar- 
ticle, which does not affect the average electron density. 
One of the difficulties in a full TCDFT calculation is that 
the electron-positron correlation functionals are poorly 
known at finite positron densities. The n+ — > limit 
used in the CONV scheme is known better. 

The thermalized positron in the lattice is in the k = 
state. When studying bulk systems we calculate the 
positron wave function at the F point (k = 0). In the case 
of positrons localized at vacancies the energy eigenvalue 
corresponding to an isolated vacancy is broadened due 
to the supercell approximation to a narrow band of en- 
ergies. We integrate over the lowest lying positron band 
by calculating the positron wave function both at the F 
point and at the Brillouin zone boundary point L and 
using the average of the respective resultsi^S 

B. Annihilation rate models 

The positron lifetime r is the inverse of the annihilation 
rate A which in a given system is proportional to the 
overlap of the electron and positron densities, 

A = i = Trr^c J dr7i_(r)n+(r)5(n_(r),7i+(r)). (2) 

Above, Te is the classical electron radius and c the speed 
of light. The enhancement factor g(n_(r), n+(r)) (the 
contact density or the electron-positron pair correlation 
function evaluated at the positron) takes into account 
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the increase in the annihilation due to the screening 
cloud of electrons around the positron. [The correspond- 
ing result obtained by omitting this factor is called the 
independent-particle model (IPM) annihilation rate.] In 
the TCDFT and within the LDA g is written as a func- 
tion of both the local electron and positron densities. In 
the CONV scheme the n+ ^ limit of the enhancement 
factor, denoted by 7(n_(r)), is used. Also Gilgien et al}^ 
used the n+ — > limit of the enhancement factor in their 
calculations but they calculated the positron density self- 
consistently within the TCDFT. This scheme has been 
shown to lead to rather localized positron states and too 
low core electron annihilation rates in comparison with 
experiments^^ 

The enhancement factor in the Boroiiski-Nieminen 
two-component formalism^ is based on the results of 
the many-body calculations by Lantto.'^" Gilgien et alm^ 



and Barbiellini et alM^ have used the 



limit 



parametrizations consistent with the correlation energy 
results of Arponen and Pajannei^ 

The LDA systematically underestimates positron life- 
times in materials because it overestimates the annihila- 
tion with core electrons for which the correlation effects 
are less important i^SiSi Therefore, Barbiellini et alm^^ 
have presented a gradient-corrected scheme in which the 
enhancement factor 7 is interpolated between the LDA 
(7 = 7LDA, zero gradient) and the IPM values (7 = 1, 
infinite gradient) as a function of the charge density gra- 
dient Vn_ . Effectively, the interpolation means that the 
annihilation with valence electrons in the interstitial re- 
gion is described using the LDA but when the density 
gradient is high (as near nuclei where the rapid oscil- 
lations of core and valence wave functions take place) 
the enhancement factor decreases and approaches the 
IPM limit (7 = 1). The interpolation form contains one 
semi-empirical parameter. The value a — 0.22 has been 
found to give with the Arponen-Pajanne enhancement 
lifetimes in good agreement with the experiment uSSsSi 
One must note that also the correlation potential for the 
positron is gradient corrected in the scheme by Barbi- 
ellini et at However, the different enhancement factors 
cause directly most of the differences compared to the 
BN-LDA. The LDA parametrization of the correlation 
potential is the same in both schemes. 



C. Schemes for the calculation of momentum 
distributions of annihilating electron-positron pairs 

The IPM formula for the momentum distribution of 
annihilating electron-positron pairs is written as 



P(P) 



J dre-*P->+(r)^j(r) 



(3) 



where "0+ (r) and tjjj (r) are wave functions of the positron 
and the electron on orbital j, respectively. The summa- 
tion goes over the occupied electron states. The IPM 



is often used because it gives, in contrast to the anni- 
hilation rate, a rather good qualitative correspondence 
(shape of the momentum distribution) with experiments. 
A common way to take into account the electron-positron 
correlation effects is to introduce in the IPM expres- 
sion the position-dependent LDA enhancement factor 
^7(n_(r)) 
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P(P) 



drt 



'"V'+(r)V'j(r)A/7(n_(r)) 



(4) 

We call this the state-independent LDA scheme. Eq. Q 
is, at least in a homogeneous system, consistent with the 
total annihilation rate A of Eq. (O . Namely, one should 
obtain A by integrating over the momentum 



A = y dpp(p). 



(5) 



The state-independent LDA scheme is motivated by the 
enhancement factor of the contact density, but it is not 
obvious how the screening really modifies the (electron) 
wave function. One can consider the position-dependent 
enhancement factor ■y/7(n_(r)) as a factor describing 
the distortion of the two-body wave function ijj^{r)ijjj{r) 
(where both the electron and the positron reside at the 
same point) due to the short-range electron-positron cor- 
relation. What is problematic is that the two-body wave 
function is distorted everywhere at the instant of the an- 
nihilation although the screening is a local phenomenon. 
This causes the correlation effects to be overestimated in 
the state-independent LDA scheme. 

In the so-called state-dependent schemeS^ a constant 
electron-state-dependent enhancement factor 7^ is used. 



P(P) = T^rlc^^jj J drt 



"■V+(r)V,(r) 



(6) 



ilPM 



The enhancement factor is written as = Xj/\j 
where Xj is the annihilation rate of the state j within 
the LDA or the GGA, 



A, 



irrlc 



dr7(n_(r))n_(. (r)nj(r), 



(7) 



and Aj-^^ is the annihilation rate within the IPM (7 = 1). 
Above, nj{v) = \^j{v)\^ is the electron density of the 
state j. In the state-dependent scheme the momentum 
density of a given (electron on a certain orbital) anni- 
hilating electron-positron pair is (apart from the factor 
7j) the same as in the IPM, i.e. the enhancement, which 
in a sense is averaged over the electron-positron pair, af- 
fects only the annihilation rate Xj not the shape of the 
momentum distribution of the orbital j. Eq. (jSJ is again 
satisfied. The problems related to the state-independent 
LDA scheme are avoided because the enhancement factor 
affects only the probability of annihilation of the positron 
with each electron state. 
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D. Projector augmented-wave method 

1. Wave functions in the projector-augmented wave method 

We use the projector augmented-wave (PAW) method'' 
to describe the valence electron wave functions in solids. 
The PAW method is a full-potential all-electron method 
related both to the pseudopotential method and to the 
linearized augmented-plane-wave method (LAPW). It is 
based on a linear transformation between all-electron 
(AE) valence wave functions |^) and soft pseudo (PS) 
valence wave functions The transformation can be 
written as (for details see Ref.0) 

m = \^)+Y.^\<t>^)-\~^^)m^). (8) 
i 

where and are AE and PS partial waves localized 
around each nucleus, and {pi\ are soft, localized projec- 
tor functions probing the local character of the PS wave 
function j^*). Index i stands for the site index i?, the an- 
gular momentum indices (?, m) and an additional index k 
referring to the reference energy £ki ■ The solution of the 
self-consistent electronic structure for a given solid sys- 
tem means the solution of the PS wave functions. They 
are represented by plane-wave expansions in the Vienna 
Ab-initio Simulation Package^^iSi^ (vASP) which we are 
using. The construction of the AE wave functions in the 
PAW method is described in detail in Ref. U The AE 
valence wave functions are orthogonal to the core states 
treated within the frozen-core approximation (free atom 
wave functions are used). 

When calculating momentum distributions of annihi- 
lating electron-positron pairs, we construct the AE wave 
functions j'l') according to Eq. © in the Fourier space 
and then Fourier transform them to the real space. In the 
case of positrons localized at defects, the summation over 
R can be limited only to the atoms surrounding the de- 
fect. The positron state is solved in the real space. Then 
the products of the positron and electron wave functions 
are calculated and Fourier transformed [see Eq. ©]. As 
a result we have a three-dimensional momentum distri- 
bution on the reciprocal lattice of the superlattice. Us- 
ing a dense k-point mesh for electron wave functions we 
decrease the lattice constant in order to increase the mo- 
mentum resolution and to get a converged result. Then, 
by integrating over the planes perpendicular to the cho- 
sen momentum distribution the Doppler spectrum is ob- 
tained with a sufficient resolution. 

It is sufficient to use a typical value of about 250- 
400 eV for the kinetic energy cutoff of the plane-wave 
expansions when calculating the PS wave functions for 
the determination of the momentum distribution. The 
momentum components of the partial waves in Eq. (jS)) 
can be taken into account up to an arbitrary value Pmax- 
We have found that the value Pmax = 70 x 10"'^ mgc is 
enough to guarantee that the Doppler spectrum (projec- 
tion of p(p) on the axis) converges up to the momen- 



tum of 40 X 10~^ ttiqc, which is usually required when 
comparing results with coincidence Doppler broadening 
experiments. 

The PAW method describes also the high-momentum 
Fourier coefficients of valence wave functions accurately, 
which is important when one compares theoretical results 
with experimental coincidence Doppler spectra. The effi- 
ciency and the flexibility of the method are also great 
benefits in the study of defects in solids. It also en- 
ables one to treat first-row elements, transition metals 
and rare-earth elements. 



2. Constructing the effective potential for the positron 

Although the PAW method is an AE method, we do 
not in practice construct the AE valence charge density n 
in the three-dimensional real-space grid when construct- 
ing the effective potential for the positron or calculating 
the total annihilation rate A. A sufficiently good approx- 
imation is to approximate n with n -I- n, where h is the 
PS valence charge density, calculated from the PS wave 
functions l^*), and n denotes the compensation charges 
as defined in Ref. 0. (Here we adopt the notation of 
Ref. 13) The compensation charges h guarantee that the 
approximate Hartree potential due to the valence elec- 
trons, vii[n + n\, is equal to the AE Hartree potential 
vii[n\ everywhere except near the nuclei, inside the local- 
ized compensation charges h [r < r^ompj where r^omp's 
are the cutoff radii of the compensation charges). The 
charge density h + fi itself is correct outside the radii 
(> '"comp)j the cutoff radu for the partial waves \(f)i) 
and from nuclei. Typically the radii are, depending 
on the element, of the order oi r\. — 1.2 ... 2.3 and 
r^Qjjj == 0.8... 2.0 ao, where oq is the Bohr radius (see 
Ref.lal. 

Our approximation is justified by the fact that near the 
nuclei the positron density is vanishingly small because 
of the Coulomb repulsion of the nuclei. Thereby, the 
positron state is not considerably affected and the overlap 
of the electron and positron densities (the annihilation 
rate A) does not appreciably change. Note, however, that 
we calculate rijir), the charge density of the state j in 
Eq. |(7J) represented in the three-dimensional real-space 
grid, directly from the AE wave function 

After calculating the Coulomb potential due to both 
the valence and core electrons and nuclei, VY^[n-\-fl-\-nzc\^ 
we calculate the electron-positron correlation potential 
and solve the positron state '0+ (r) in the effective poten- 
tial V:,_(r) of Eq. (P) using the RQMG solver This is a 
fast and accurate method for our purpose where only the 
positron ground state corresponding to a rather smooth 
wave function in the interstitial region has to be solved. 
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E. Description of positron annihilation with core 
electrons 



When modeling the positron annihilation and calcu- 
lating momentum distributions of annihilating electron- 
positron pairs we describe the core electrons and the core 
electron charge density using atomic orbitals of isolated 
atoms calculated within the DFT and the LDA. In the 
calculation of the electron-positron pair wave function for 
the core electron Dopplcr spectrum we use an isotropic 
parametrized positron wave function* of the form 



iP+{r)^C{ai + [eTi{r/a2)r}, 



(9) 



where C is a normalization factor and ai,a2,a3 are pa- 
rameters determined by fitting Eq. to a spherically- 
symmetric positron wave function calculated with the 
LMTO method within the atomic-spheres approxima- 
tion (ASA). It is sometimes questionable to assume 
the positron wave function to be spherically symmet- 
ric around the nuclei when calculating the core elec- 
tron Doppler spectrum. In the perfect bulk the positron 
wave function is very isotropic close to the nuclei. For 
the positron trapped by a vacancy-type point defect 
the decay of the positron wave function in the neigh- 
boring ion cores is similar to that in the bulk. The 
anisotropy around the nuclei causes extra localization in 
the positron-core-electron overlap, which causes some in- 
crease of the positron momentum which is omitted in 
the model. However, this is expected to be small in com- 
parison with the positron momentum due to the decay 
toward the nuclei. 

To test the effects of the frozen-core approximation 
and the isotropic positron wavefunction used when mod- 
eling the core-electron Doppler spectrum we have made 
two calculations for As vacancy (Vas) in GaAs by treat- 
ing the 3d electrons of Ga first as valence electrons and 
then as core electrons. In the former calculation the 
full three-dimensional positron wave function is used in 
constructing the three-dimensional electron-positron pair 
wave functions corresponding to the Ga 3d electrons 
whereas in the latter calculation the positron wave func- 
tion has the isotropic form of Eq. I^. The intensities of 
the results differ at high momenta due to the different de- 
gree of self-consistency but when one compares VAs/bulk 
ratios of Doppler spectra the results coincide. (As long 
as the 3d electrons of Ga are treated consistently in the 
bulk and defect calculations.) 



F. Calculation of forces on ions due to a localized 
positron 

When modeling relaxations of ions around vacancy de- 
fects the effects of the localized positron are included 
by using the so-called atomic-superposition (ATSUP) 
approximation^ in which the charge density and the 
Coulomb potential are constructed from those of isolated 
atoms. In the CONV scheme the total energy is the sum 



of the total energy of the electron-ion system and the 
positron energy eigenvalue e+. Then the force due to the 
positron on ion j is the negative gradient of the positron 
energy eigenvalue e+ with respect to the position of the 
ion Rj. According to the Hellman-Feynman theorem, 

F+ = -V,e+ = -Vj(V'+|i?|?A+) = -(?A+|V,i/|V+), 

(10) 

where the positron wave function is assumed to 

be properly normalized. Within the ATSUP approxima- 
tion and the LDA the effective potential for the positron 
[Eq. (P)] is of the form 



V+{r) 



E 



R, 



|) + Korr( ^"-''(1^" 



Ri 



(11) 

where Vq^'^^ and '•' are the Coulomb potential and the 
charge density of the free atom j, respectively. Inserting 
this into Eq. H1U|) gives for the force 



dr n+ (r) 



dn^''' (r) 



dn 



dr 



dr 



|r-R,| 



(12) 



|r-R,| 

r-R, 



R.I 



The calculation of the positron-induced forces is fast 
within the ATSUP method. The forces are used with 
calculated electron-ion and ion-ion forces in order to find 
the relaxed ionic configuration of the defect. The ATSUP 
method itself does not give the possibility to study the 
relaxation of charged defects. However, if the positron 
state is solved in the self-consistent Coulomb potential 
instead of the superimposed potential of Eq. lll|l . the 
effect of the charge state comes into play, for example, 
in the case of negatively charged vacancies in semicon- 
ductors, via the stronger localization of the positron in 
comparison with the neutral vacancy, which will result 
in larger positron-induced forces. The approximations 
made when using Eq. H12() are tested below and com- 
pared with TCDFT results. 



III. PERFECT BULK SYSTEMS 

A. Testing the PAW method 

We begin the testing of the PAW method for the 
calculation of momentum distributions of annihilating 
electron-positron pairs by comparing the results to those 
of the ATSUP approximation where atomic orbitals are 
used also for valence electron states.^ Here we also 
demonstrate the effect of the PAW transformation on the 
PS wave functions by showing also results calculated us- 
ing only the PS wave functions of the PAW method. In 
the calculation of the electronic structures we employ the 
LDA exchange-correlation energy. The positron states 
and the annihilation characteristics are calculated within 
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the BN-LDA.^ The Doppler spectra are calculated using 
the state-dependent schemed The momentum distribu- 
tions are finally convoluted with a Gaussian function with 
a full width at half maximum (FWHM) corresponding to 
the resolution of the Doppler experiment. 

We have found out that a good way to compare theo- 
retical and experimental Doppler spectra is to plot ratios 
to a reference spectrum. This way most of the systematic 
error is canceled. For the theoretical spectra this means 
especially that the overestimation of the core annihila- 
tion in the LDA is not a problem. In Fig. ^ we show our 
results for the ratios Cu/Al and Ag/Al calculated with 
the AE-PAW method, the PS wave functions of the PAW 
method and the ATSUP method compared to the exper- 
imental data from Ref. In these calculations the Cu 
3d and Ag 4d electrons are treated as valence electrons 
in the PAW calculations. 

Due to the non-selfconsistent construction of the va- 
lence electron wave functions the ATSUP method gives 
good results only at high momenta where the annihila- 
tion with the core electrons dominates. The results ob- 
tained with the AE-PAW method are clearly better in 
the low-momentum region of the spectra. In the high- 
momentum region these two results are equally good. 
The better compatibility of the crude ATSUP approxi- 
mation with the experiment at high momenta in Fig.^a) 
may be just a coincidence. The PS wave functions of the 
PAW method fail to represent the high-frequency oscil- 
lations of the valence wave functions, especially those of 
the Cu 3d and Ag 4d electrons, in the core region al- 
though they predict the ratios well up to the momentum 
of about 10 X 10-3 niQC. The AE-PAW resuhs in Fig. □ 
are practically combinations of the ATSUP (at high mo- 
menta) and the PS-PAW (at low momenta) results. Be- 
cause the quality of the PS-PAW results is comparable 
to calculations made using norm-conserving pseudowave 
fimctions^''-^^ the tests made in this section clearly show 
the benefits of the use of the PAW method in the accu- 
rate calculation of momentum distributions of annihilat- 
ing electron-positron pairs. We also note that positron 
lifetimes obtained for different bulk systems agree per- 
fectly with previous all-electron results calculated with 
the same enhancement factor. 



B. State-dependent scheme vs. state-independent 
LDA scheme 

In this section we compare the two above-mentioned 
ways to take into account the electron-positron correla- 
tion in the calculation of the momentum distribution of 
annihilating electron-positron pairs: the state-dependent 
scheme of Eq. © and the state-independent LDA scheme 
of Eq. For simplicity, we use only bulk materials (Cu, 
Ag, Al, Si, Mo, and Fe) as examples. We use the experi- 
mental lattice constants and for Fe the experimental bcc 
structure. From now on we use the AE-PAW method 
within the frozen-core approximation. 
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FIG. 1: (Color online) (a) Bulk Cu / bulk Al and (b) bulk 
Ag / bulk Al ratio curves of momentum distributions of an- 
nihilating electron-positron pairs calculated using the state- 
dependent scheme. The experimental data— is shown with 
circles. The theoretical curves are convoluted with a Gaus- 
sian function with a FWHM of 4.3 x 10"^ moc. 
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FIG. 2: (Color online) Bulk Cu / bulk Al ratio curve of mo- 
mentum distributions of annihilating electron-positron pairs. 
The comparison between the state-dependent and the state- 
independent LDA schemes is shown. The experimental data 
is from Ref. [33 . The theoretical curves are convoluted with a 
Gaussian function with a FWHM of 4.3 x 10"^ moc. 



We show in Figs. [21 and the Doppler spectra of Cu, 
Ag and Al calculated within the both above-mentioned 
schemes in the logarithmic scale and also normalize 
the Doppler spectra to the one of Al. The theoretical 
spectra and ratios are compared with the experimental 
ones by Nagai tt alm^ Both schemes describe the low- 
momentum region due to the annihilation with valence 
electrons well but the ratios tell, as seen before,* that the 
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FIG. 3: (Color online) Bulk Ag / bulk Al ratio curve of mo- 
mentum distributions of annihilating electron-positron pairs. 
The comparison between the state-dependent and the state- 
independent LDA schemes is shown. The experimental data 
is from Ref. Issi . The theoretical curves are convoluted with a 
Gaussian function with a FWHM of 4.3 x 10~^ moc. 

state- independent LDA scheme fails to describe the high- 
momentum region of the Doppler spectra which arises 
from the annihilation with core electrons. The ratio 
Cu/Al (see Fig.|21l calculated using the state-independent 
LDA scheme is not even in a qualitative agreement with 
the experiment at high momenta. On the contrary, the 
state-dependent scheme describes the ratios quite accu- 
rately in both Figs. [3 and 13 However, when one looks at 
the absolute values of the spectra, the intensities at high 
momenta are in better agreement with the experiment in 
the results calculated using the state-independent LDA 
scheme but there is some unphysical oscillation in the 
spectra that does not exist in the state-dependent LDA 
results. 

Further examples are shown in Figs. ^ and |S1 where 
the Doppler spectra of Al, Si, Mo and Cu are normal- 
ized to the one for Fe. For Fe we consider its magnetic 
ground state. The state-independent LDA fails again at 
high momenta, the result for Cu being again in the worst 
agreement with the experiment. 



C. LDA vs. GGA within the state-dependent 
scheme 

In Figs. 01 13 El and[7|we show similar comparisons be- 
tween the BN-LDA and the GGA by Barbiellini et al^^ 
The state-dependent scheme is used. The results calcu- 
lated using the BN-LDA for the annihilation rates are in 
a better agreement with the experiment. In Figs.Eland[7| 
the GGA tends to give too high values for the ratio at 
high momenta. In these particular examples the failure 
is mainly due to the large decrease of the relative core an- 
nihilation rate of Al calculated with the GGA compared 
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FIG. 4; (Color online) Ratio curves of momentum distri- 
butions of annihilating electron-positron pairs calculated us- 
ing different approximations. The experimental data is from 
Ref. |3^ The theoretical curves are convoluted with a Gaus- 
sian function with a FWHM of 4.7 x 10~^ moc. 




p^(10^m„c) 

FIG. 5: (Color online) Bulk Cu / bulk Fe ratio curve of mo- 
mentum distributions of annihilating electron-positron pairs 
calculated using different approximations. The experimental 
data is from Ref. Is^ . The theoretical curves are convoluted 
with a Gaussian function with a FWHM of 4.7 x 10~^ moc. 



to the one according to the BN-LDA. Because the state- 
dependent scheme is used the shapes of the contribu- 
tions due to individual orbitals are the same and also the 
shapes of the ratios are very similar. The ratios Al/Fe, 
Mo/Fe and Cu/Fe calculated with the GGA are shown 
in Figs. 0] and O The results are in a rather good agree- 
ment with the experiment and only slightly worse than 
the ones calculated with the BN-LDA state-dependent 
scheme. The relative underestimation of the core annihi- 
lation in Al can be seen also in the ratio Al/Fe, which is 
lower at high momenta than the BN-LDA result. In con- 
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FIG. 6: (Color online) Bulk Cu / bulk Al ratio curve of mo- 
mentum distributions of annihilating electron-positron pairs. 
The comparison between the BN-LDA and the GGA for an- 
nihilation rates within the state-dependent scheme is shown. 
The experimental data is from Ref. 35» The theoretical curves 
are convoluted with a Gaussian function with a FWHM of 
4.3 X 10"=* moc. 
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FIG. 7: (Color online) Bulk Ag / bulk Al ratio curve of mo- 
mentum distributions of annihilating electron-positron pairs. 
The comparison between the BN-LDA and the GGA for an- 
nihilation rates within the state-dependent scheme is shown. 
The experimental data is from Ref.lSa. The theoretical curves 
are convoluted with a Gaussian function with a FWHM of 
4.3 X 10"^ moc. 



trast, as seen in Figs. Eland 13 the GGA describes better 
the absolute intensities of the Doppler spectra because 
the annihilation rates of core orbitals are decreased. 

The failure of the GGA in the ratios of Doppler spec- 
tra can be traced back to the semiempirical interpola- 
tion form of the GGA enhancement factor. Although its 
zero- and high-gradient limits are well defined the inter- 
polation form is only an approximation. Moreover, the 



free parameter a is just fixed to give lifetimes that are in 
good agreement with the experiment. Clearly, the BN- 
LDA succeeds to describe better the relative magnitudes 
of the annihilation rates Xj between different electronic 
states and different elements. 

We conclude that our choice for the approximation to 
be used with the state-dependent scheme is the BN-LDA 
because it gives better results than the GGA when com- 
paring intensity ratios with the experiment. It is also 
simpler and more justifiable. 



IV. VACANCY DEFECTS IN SOLIDS 

The following step is to demonstrate that our scheme 
works also for positrons localized at vacancy-type de- 
fects. The CONY scheme has been shown to yield for 
a given ionic structure lifetimes^^ and also other anni- 
hilation characteristics^^ in good agreement with two- 
component calculations based on the Borohski-Nieminen 
formalism. 

In this section we compare our results for vacancies 
in metals and semiconductors to experimental Doppler 
broadening results. We also study the effect of the 
positron- induced forces on ions neighboring vacancies. 
Further, we compare the ionic relaxations to previous 
two-component results and compare Doppler spectra cal- 
culated with the relaxed structures with experimental 
ones. 



A. Relaxation tests 

We study the effect of the localized positron on the 
relaxation of ions surrounding a vacancy by calculating 
the ionic structures of monovacancies in bulk Si, Al and 
Cu with and without the localized positron. We consider 
first only isotropic relaxations. For Si we use a cubic 64- 
atom supercell and for Al and Cu cubic 108-atom super- 
cells. In the electronic structure calculations we sample 
the Brillouin zone using 4'^, 8'^ and 6^ Monkhorst-Pack^ 
k-point meshes, respectively. Self-consistent LDA lattice 
constants are used in all calculations. We use the value 
0.01 eV/A as a stopping criterion for the forces on ions 
when finding the self-consistent ionic configurations. 

We consider two different approximations for the force 
calculation. We solve the positron state either in the 
PAW potential or in the ATSUP potential. The first 
approximation is better in the sense that the positron 
density is able to follow the changes in the electronic 
structure but our force expression, Eq. (|12|l . is based on 
the ATSUP approximation and therefore the PAW po- 
tential is not consistent with the potential of the force 
calculation and thus the total energy minimum does not 
correspond to vanishing forces. In contrast, the latter 
more crude approximation is consistent with the force 
expression used and gives thus a well-defined and consis- 
tent total energy minimum. 
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TABLE I: Relaxations and corresponding positron liletimes 
T for monovacancies in different bulk materials. A positive 
(negative) number denotes the isotropic outward (inward) re- 
laxation in percentage of the nearest neighbor distance with 
respect to the unrelaxed (ideal) vacancy. The table includes 
results calculated without the effect of the positron and with 
the positron state solved in the PAW/ATSUP potential. The 
relaxation is restricted to the symmetric breathing-mode re- 
laxation. The numbers in parenthesis are calculated using 
a larger 216-atom supercell. Computed bulk lifetimes are 
208 ps, 159 ps and 95 ps for Si, AI and Cu, respectively. 





no 


e+ 


PAW 


ATSUP 




rel. (%) 


r (ps) 


rel. (%) r (ps) 


rel. (%) r (ps) 


Si 


-10.4 


215 


4-5.6 (-1-11.3) 256 (272) 


+5.9 257 


Al 


-1.7 


219 


+2.8 242 


+2.9 251 


Cu 


-1.3 


146 


+2.4 163 





The relaxations obtained with and without the 
positron are listed in Table U The effect of the positron 
on the relaxations is clear; in all of these examples the in- 
ward relaxation is transformed to an outward relaxation 
due to the positron. The relaxations obtained with the 
PAW and the ATSUP potentials are very similar. We 
have studied the Si vacancy with the localized positron 
using also a larger cubic 216-atom supercell and only the 
r point. The results obtained are shown in Table Q in 
parenthesis. The larger outward relaxation is explained 
by the fact that in the larger supercell the ions can relax 
more freely; interactions between periodic images of the 
vacancy are not as dominant as in the 64-atom super- 
cell. In Fig. |S1 we have plotted different energy compo- 
nents of the Vsi and Vai systems as functions of vacancy 
relaxation relative to the relaxed geometry (ionic struc- 
ture obtained including the effect of the forces due to the 
positron). Note, that one can not compare absolute ener- 
gies between the results calculated with different poten- 
tials for the positron. The energy minimum for the PAW 
positron potential plus potential due to the electron-ion 
system is in both systems at about 1 % smaller relaxation 
than the structure given by the relaxation (zero forces). 

Our computational positron lifetimes for the monova- 
cancies in Al and Cu are in good agreement with the 
experiment. The exp erimental lifetimes for Vai and Vcu 
are 251 ps (Ref. l38l) and 179 ps (Ref. Is^ . respectively. 
They are high in comparison with the computed lifetimes 
since we use the LDA enhancement factor but we find a 
good agreement when we compare the vacancy-bulk life- 
time differences with the experiment. (The experimental 
lifetimes for bulk Al and bulk Cu are 170 ps and 120 ps, 
respectively.^^) The calculated and experimental results 
for the monovacancy in Si are compared in section FlV CI 




Relaxation (%) 

FIG. 8: (Color online) Different energy components as func- 
tions of vacancy relaxation in % of the nearest-neighbor dis- 
tance (relative to the relaxed geometry, only nearest neighbor 
ions are moved) for Vsi in Si and Vai in AI. The solid lines 
denote the total energy, the dashed lines the energy of the 
electron-ion system, and the dash-dotted line the positron en- 
ergy eigenvalue. The figures show results in which the forces 
are calculated using a positron state solved in the PAW po- 
tential (o) and in the ATSUP potential (o). 



B. Ga vacancy in GaAs 

The triply negative Ga vacancy in GaAs has been ex- 
tensively studied by Puska et alw^ using the TCDFT and 
different schemes (including the CONV) for the electron- 
positron correlation. Furthermore, also experimental co- 
incidence Doppler broadening data existsr^- Thus, us- 
ing the Ga vacancy as a benchmark system, it is possi- 
ble to compare simultaneously the relaxations and life- 
times obtained to two-component results and lifetimes 
and Doppler spectra to the experiment. We model the 
GaAs lattice using a cubic 64-atom supercell and sample 
the Brillouin zone in the electronic-structure calculations 
with a 4"^ Monkhorst-Pack k-point meshi^ In the case of 
the triply negative charge state of Vca all the localized 
states in the band gap are occupied and there is no sym- 
metry lowering Jahn- Teller relaxation. By including the 
effect of the positron-induced forces on the relaxation we 
get an inward relaxation of 5.9 % in very good agreement 
with the inward relaxation of 6.6 % previously obtained 
using the two-component BN schemeii^ In the work by 
Puska et al^ the authors did not calculate relaxations 
using the CONV scheme but showed that the total energy 
curves calculated as a function of the breathing-mode re- 
laxations of the ions neighboring the vacancy using the 
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FIG. 9; (Color online) Theoretical ratio curves for the Ga 
vacancy in GaAs. The experimental data, measured from 
electron-irradiated GaAs, is from Ref. l40l The theoretical 
curves are convoluted with a Gaussian function with a FWHM 
of 5.5 X 10~^ moc. 



TABLE II: Relaxations and lifetimes for different charge 
states of Vca in GaAs. The numbers in parenthesis are calcu- 
lated using a larger 216-atom supercell. The relaxation is re- 
stricted to the symmetric breathing-mode relaxation. A pos- 
itive (negative) number denotes isotropic outward (inward) 
relaxation. The computed bulk GaAs lifetime is 208 ps. 



defect rel. (%) r (ps) r — rbuik (ps) 

unrel. 0.0 249 38 

Vq-; -2.8 (-2.3) 237 (244) 26 (33) 

V^; -5.9 229 18 

Vca (exp^) 260 30 



the two-component BN and the CONV schemes nearly 
coincide. When this and the good agreement in relax- 
ations are taken into account, one can conclude that our 
scheme for the calculation of forces gives very similar 
results to ones calculated using the two-component BN 
formalism. 

In Fig. 1^ we show the computed Doppler spectra (nor- 
malized to that of bulk GaAs) obtained for the relaxed 
structures of Vca in the charge states 3— and 1— com- 
pared with the experiment and the computed one for the 
neutral, ideal (unrelaxed) Vca- The corresponding re- 
laxations and lifetimes are tabulated in Table ^ Only 
isotropic relaxations have been considered using the small 
64-atom supercell although a symmetry-breaking relax- 
ation is expected for the 1— charge state. The agreement 
with the experiment both in the Doppler spectrum and 
in the lifetime (relative to the bulk one) is best for the 
1— charge state. The inward relaxation for the 3— is too 
strong compared with the experiment because the life- 
time and the ratio curve in Fig. |5| are too close to the 
bulk values. 

We have also calculated the relaxation of Vq^ using 
a larger cubic 216-atom supercell (with this supercell 
we use only the F point and do not treat the Ga 3d 



electrons as valence electrons). The calculation gives 
slightly smaller inward relaxation than that with the 64- 
atom cell. The results are given in Table ITU in parenthe- 
sis. We also break the Td symmetry by displacing the 
nearest-neighbor atoms of the vacancy in order to create 
a symmetry-breaking relaxation with the expected C21, 
symmetry. The corresponding lifetime is 243 ps, which 
is almost the same as when assuming the T^ symmetry. 
(Note that Fig. |21 does not include any results calculated 
with the 216-atom supercell.) 

The relative core annihilation rate (the experimentally 
measured relative W parameter, relative wrt bulk value) 
is sensitive to the treatment of the electron-positron cor- 
relation effects. The correlation potential used affects 
the degree of the localization of the positron and thereby 
the positron-core electron overlap and the core annihila- 
tion rate.^^ Puska et ali^ obtained for the triply negative 
Ga vacancy in GaAs the values of 0.88 and 0.34 using the 
two-component BN formalism and the scheme by Gilgien 
et al.^ respectively. The computational values in the 
present work (estimated from Fig. O are of the order of 
0.8 in good agreement with the experiment and the previ- 
ous result obtained with the BN formalism. Our scheme 
has already previously been successful in describing the 
relative core annihilation rates in the case of Ga vacancy 
in GaN^^ and monovacancy in Ali24 



C. Neutral monovacancy in Si 

The neutral monovacancy in Si is a very complex sys- 
tem with a flat potential-energy surface and several com- 
peting local minima as a function of ion positions 
One can get even qualitatively different results with two 
different approximations e.g. for the exchange-correlation 
potential)^ We study now the monovacancy in Si in- 
cluding the forces caused by the localized positron. We 
use the 216-atom supercell and the F point, begin from 
the structure obtained as a result from the isotropic re- 
laxation (see Table ^ and break the Td symmetry by 
displacing the nearest-neighbor atoms of the vacancy in 
order to create a symmetry-breaking relaxation with the 
expected D2d symmetry. We confirm the fact pointed 
out by several earlier studies^S^Siilii^ that a supercell 
with at least 216 atoms is needed in order to obtain a 
convergence with respect to the supercell size. 

Our calculation gives an outward-relaxed structure 
with a slight D2d symmetry. The corresponding lifetime 
is 270 ps; only 2 ps less than that for the vacancy con- 
strained to the Td symmetry. In Fig. 1101 we show the 
obtained vacancy/bulk ratio of Doppler spectra. No ex- 
perimental Doppler broadening data exists for the mono- 
vacancy in Si. Makinen et aZ4& and Polity et al.*'' ob- 
tained the positron lifetimes of 273 and 282 ps, respec- 
tively, for the monovacancy in Si created by electron ir- 
radiation. Taking into account that the calculated bulk 
lifetime with the LDA lattice constant used (208 ps) is 
lower than the experimental ones by Makinen et al. and 
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FIG. 10: Theoretical ratio curve for the neutral Si vacancy in 
Si. The the data is convoluted with a Gaussian function with 
a FWHM of 3.7 x 10"^ moc 

Polity et al. (221 and 218 ps, respectively) we conclude 
that our result for the lifetime is in a good agreement 
with experiment. 

Previously, Saito and Oshiyam£4£ and Makhov and 
Lewis^ have studied vacancies in Si within the two- 
component scheme by Gilgien et ali^ including the ef- 
fects of forces due to the positron. Their lifetimes for the 
monovacancy are reasonable but the relative W param- 
eter of 0.28 estimated from Saito and Oshiyama's data 
is very low in comparison with our result, which is 0.72 
(both evaluated using calculated annihilation rates as in 
Ref.iia). 

V. SUMMARY AND CONCLUSIONS 

In conclusion, we have presented an accurate scheme 
for the calculation of momentum distributions of annihi- 
lating electron-positron pairs in solids based on the pro- 
jector augmented- wave method. 

We have compared three commonly used approaches 
for the momentum distributions within the DPT frame- 
work. We have shown that the most appropriate way 
to take into account screening effects in the calculation 
of momentum distributions is to use a constant, state- 



dependent enhancement factor. Further, we have demon- 
strated that a position-dependent enhancement factor 
gives unphysical results when ratios of Doppler spectra 
are considered. The differences in results of the BN- 
LDA and the GGA by Barbiellini et al. is not large. 
Except for one of the studied elements (Al) the GGA 
gives comparable results. We choose to rely on the BN- 
LDA because of its simplicity and lack of semi-empiric 
parameters. We underline that our choices are based on 
the theory-experiment comparison of ratios of momen- 
tum distributions for different materials rather than their 
absolute values. The latter may often be better described 
in the GGA and in the position-dependent enhancement 
schemes. 

In addition to bulk solids our scheme is also reliable in 
the case of defects in materials. The comparison of our re- 
sults for the Ga vacancy in GaAs to experiment suggests 
that the Ga vacancy seen in the experiment is negative 
but less than triply negative, which is the charge state 
suggested by a recent ab initio study^S Por the neutral 
monovacancy in Si we have presented a prediction to be 
compared with future lifetime and coincidence Doppler 
broadening experiments. 
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